PRO MAKE_Z_PLOT
  
  ;; Define the COMMON BLOCK
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  FORWARD_FUNCTION BGPS_FFORE_MODEL, BGPS_FFORE_LOOPIE, STAR_DENSITY
  RESOLVE_ROUTINE, 'BGPS_FFORE', /COMPILE_FULL_FILE, /EITHER, /NO_RECOMPILE
  
  
  IF n_elements(rb3) EQ 0 THEN restore,'./irdc_dist_model/bgps_rb3.sav',/ver
  restore,'./irdc_dist_model/data/case2_irdc_TdFIT.sav',/ver
  n = n_elements(rb3)
  
  omni_lbd2rz,rb3.l,rb3.b,fitdist,R,Z
  
  ind = where(constrain, nind) & print,nind
  rb3 = rb3[ind]
  R = R[ind]
  Z = Z[ind]
  n = nind
  
  ;;===========================================
  ;; Make the PS plot
  
  myps,'./irdc_dist_model/analysis_plots/z_position_subs.eps'
  
  ;; DIVIDE the group into 3 sections by Rgal
  rs = r[sort(r)]/1.d3
  zs = z[sort(r)]
  
  cut = [0,long((findgen(3)+1.)/3.*n)-1]
  print,rs[cut],format="(F0.2)"
  
  pm = ' '+cgSymbol('+-')+' '
  xrange = [-150,150]
  
  ;; hist = HISTOGRAM(z, binsize=5, loc=zout)
  plothist,zs[cut[0]:cut[1]],zout,hist,psym=10,ytit='N per 5 pc bin',thick=3,$
           xtit='Z (Vertical Position) [pc]',xr=xrange,/xst,bin=5,charsize=1.0
  yfit = mpfitpeak(zout, hist, A, NTERMS=3, ERROR=sqrt(float(hist)),PERROR=perr)
  oplot,zout,yfit,color=cgColor('Black'),thick=3
  txtstr = ['Fit centroid = '+string(A[1],pm,perr[1],format=$
                                     "(F+0.1,A0,F0.1)")+' pc',$
            'FWHM = '+string(A[2]*sqrt(8. * alog(2.)),format="(F0.1)")+' pc',$
            'Gauss Sig = '+string(A[2],pm,perr[2],format="(F0.1,A0,F0.1)")+' pc',$
            'R!dgal!n = '+string(rs[cut[0]],format="(F0.2)")+' - '+$
            string(rs[cut[1]],format="(F0.2,' kpc')")]
  
  al_legend,txtstr,box=0,/top,/left,charsize=1.0,spacing=1.5
  
  
  plothist,zs[cut[1]+1:cut[2]],zout,hist,psym=10,bin=5,charsize=1.0,/overplot,$
           color=cgColor('DeepPink'),thick=3
  yfit = mpfitpeak(zout, hist, A, NTERMS=3, ERROR=sqrt(float(hist)),PERROR=perr)
  oplot,zout,yfit,color=cgColor('DeepPink'),thick=3
  txtstr = ['Fit centroid = '+string(A[1],pm,perr[1],format=$
                                     "(F+0.1,A0,F0.1)")+' pc',$
            'FWHM = '+string(A[2]*sqrt(8. * alog(2.)),format="(F0.1)")+' pc',$
            'Gauss Sig = '+string(A[2],pm,perr[2],format="(F0.1,A0,F0.1)")+' pc',$
            'R!dgal!n = '+string(rs[cut[1]+1],format="(F0.2)")+' - '+$
            string(rs[cut[2]],format="(F0.2,' kpc')")]
  
  al_legend,txtstr,box=0,/top,/right,charsize=1.0,spacing=1.5,textcolor='DeepPink'
  
  
  plothist,zs[cut[2]+1:cut[3]],zout,hist,psym=10,bin=5,charsize=1.0,/overplot,$
           color=cgColor('DodgerBlue'),thick=3
  yfit = mpfitpeak(zout, hist, A, NTERMS=3, ERROR=sqrt(float(hist)),PERROR=perr)
  oplot,zout,yfit,color=cgColor('DodgerBlue'),thick=3
  txtstr = ['Fit centroid = '+string(A[1],pm,perr[1],format=$
                                     "(F+0.1,A0,F0.1)")+' pc',$
            'FWHM = '+string(A[2]*sqrt(8. * alog(2.)),format="(F0.1)")+' pc',$
            'Gauss Sig = '+string(A[2],pm,perr[2],format="(F0.1,A0,F0.1)")+' pc',$
            'R!dgal!n = '+string(rs[cut[2]+1],format="(F0.2)")+' - '+$
            string(rs[cut[3]],format="(F0.2,' kpc')")]
  
  al_legend,txtstr,box=0,/bottom,/left,charsize=1.0,spacing=1.5,$
            textcolor='DodgerBlue'
  
  axis,xaxis=0,xthick=3,xr=xrange,/xst,xtickformat='blank_axis'
  
  myps,/done
  
  
  ;;===========================================
  ;; Make the PS plot
  
  myps,'./irdc_dist_model/analysis_plots/z_position.eps',xsize=9.5,/cmyk
  
  
  multiplot,[2,1],/doyaxis
  xrange=[-149.5,149.5]
  yrange=[0,60]
  
  plothist,z,zout,hist,psym=10,xtit='Z (Vertical Position) [pc]',$
           ytit='N per 5 pc bin',thick=3,xr=xrange,/xst,bin=5,charsize=1.0,$
           title='Full Sample',yr=yrange,/yst,ytickformat="(I0)"
  yfit = mpfitpeak(zout, hist, A, NTERMS=3, ERROR=sqrt(float(hist)),PERROR=perr)
  
  zc = findgen(1001)/5. - 100
  yc = gauss1(zc,[A[1],A[2],A[0]],/peak)
  oplot,zc,yc,color=cgColor('BLK4'),thick=6
  plothist,z,zout,hist,psym=10,xtit='Z (Vertical Position) [pc]',$
           ytit='N per 5 pc bin',thick=3,xr=xrange,/xst,bin=5,charsize=1.0,$
           title='Full Sample',yr=yrange,/yst,ytickformat="(I0)",/overplot
  txtstr = ['Fit centroid = '+string(A[1],pm,perr[1],format=$
                                     "(F+0.1,A0,F0.1)")+' pc',$
            'FWHM = '+string(A[2]*sqrt(8. * alog(2.)),format="(F0.1)")+' pc',$
            'Gauss Sig = '+string(A[2],pm,perr[2],format="(F0.1,A0,F0.1)")+' pc',$
            'R!dgal!n = '+string(rs[cut[0]],format="(F0.2)")+' - '+$
            string(rs[cut[3]],format="(F0.2,' kpc')")]
  
  al_legend,txtstr,box=0,/top,/left,charsize=0.8,spacing=1.3
  al_legend,['N = '+string(n_elements(rb3),format="(I0)")],box=0,/top,/right,$
            charsize=0.8
  
  axis,xaxis=0,xthick=3,xr=xrange,/xst,xtickformat='blank_axis'
  
  multiplot
  
  
  ;;===========================================
  ;; Make the panel for GRS-YES sources only!
  
  
  ind = WHERE(rb3.grs EQ 1 AND rb3.grs_match EQ 1, nind)
  print,nind
  
  
  
  
  plothist,z[ind],zout,hist,psym=10,xtit='Z (Vertical Position) [pc]',$
           ;;ytit='N per 5 pc bin',$
           thick=3,xr=xrange,/xst,bin=5,charsize=1.0,$
           title='GRS Distance Matched',yr=yrange,/yst
  yfit = mpfitpeak(zout, hist, A, NTERMS=3, ERROR=sqrt(float(hist)),PERROR=perr)
  yc = gauss1(zc,[A[1],A[2],A[0]],/peak)
  oplot,zc,yc,color=cgColor('BLK4'),thick=6
  txtstr = ['Fit centroid = '+string(A[1],pm,perr[1],format=$
                                     "(F+0.1,A0,F0.1)")+' pc',$
            'FWHM = '+string(A[2]*sqrt(8. * alog(2.)),format="(F0.1)")+' pc',$
            'Gauss Sig = '+string(A[2],pm,perr[2],format="(F0.1,A0,F0.1)")+' pc',$
            'R!dgal!n = '+string(min(r[ind]/1.d3),format="(F0.2)")+' - '+$
            string(max(r[ind]/1.d3),format="(F0.2,' kpc')")]
  plothist,z[ind],zout,hist,psym=10,xtit='Z (Vertical Position) [pc]',$
           ;;ytit='N per 5 pc bin',$
           thick=3,xr=xrange,/xst,bin=5,charsize=1.0,$
           title='GRS Distance Matched',yr=yrange,/yst,/overplot
  
  al_legend,txtstr,box=0,/top,/left,charsize=0.8,spacing=1.3
  al_legend,['N = '+string(nind,format="(I0)")],box=0,/top,/right,charsize=0.8
  
  axis,xaxis=0,xthick=3,xr=xrange,/xst,xtickformat='blank_axis'
  
  myps,/done,/mp
    
END
